      SUBROUTINE CNTLNS (X,Y,Z,N,ISMOPT,IEXP,JEXP,NCNTRS,CLIST,
     *                   EPSLON,IERR)
C
C     ------------------------------------------------------------------
C
C     DRIVER PROGRAM FOR COMPUTING AND DRAWING CONTOUR LINES OF
C     CONSTANT Z FOR THE FUNCTION Z = F(X,Y).
C
C
C     ARGUMENT LIST DEFINITIONS -
C
C        X       = INPUT LIST OF X VALUES
C        Y       = INPUT LIST OF Y VALUES
C        Z       = INPUT LIST OF Z VALUES
C        N       = INPUT SPECIYING THE NUMBER OF VALUES IN X,Y AND Z
C        ISMOPT  = SMOOTHING OPTION FLAG (0=NO/OFF, 1=YES/ON)
C        IEXP    = I EXPONENT VALUE FOR SMOOTHING
C        JEXP    = J EXPONENT VALUE FOR SMOOTHING
C        NCNTRS  = NUMBER OF CONTOUR LINES TO BE DRAWN
C                  (SELF COMPUTING IF NCNTRS.LE.0)
C        CLIST   = LIST OF CONSTANT CONTOUR VALUES IF NCNTRS.GT.0
C        EPSLON  = ERROR FUNCTION (NORMALIZED VALUE) RETURNED TO
C                  CALLER IF ISMOPT IS NON-ZERO
C        IERR    = RETURN ERROR FLAG
C                  = 0 FOR NORMAL RETURN
C                  = 1 FOR INVALID VALUE FOR N
C                  = 2 FOR NUMBER OF ISMOPT COEFFICIENTS GREATER
C                      THAN 'MAXCOF' OR N
C
C        (NOTE / IF NCNTRS.LE.0, THEN CLIST(1) = BASE VALUE,
C         AND CLIST(2) = INCREMENT VALUE (DELTA) )
C
C
C     ------------------------------------------------------------------
C
C
      DIMENSION  X(N),Y(N),Z(N),CLIST(2)
      DIMENSION  ZNEW(500)
      DIMENSION  IE(1494,2),ITE(1494,4),XI(1494),ETA(1494),
     *           LAMBDA(1494),IBE(1494)
      DIMENSION  IPOWR(23),JPOWR(23),COEF(23)
C
      DATA       MAXCOF /23/
      DATA       MAXPTS /500/
C
C
C        (A)
C        INITIALIZE LOCAL VARIABLES
C        AND CHECK INPUTS FOR ERRORS
C
      IERR = 0
      EPSLON = 0.0
      IF (N.LT.3.OR.N.GT.MAXPTS) GOTO 997
C
C
C        (B)
C        CALL SUBROUTINE TRAING TO TRIANGULATE X-Y DATA POINTS
C
      CALL TRIANG (X,Y,N,LEDGES,IE,IBE,ITE)
C
C        (C)
C        SMOOTHING REQUIRED? . .
C
      IF (ISMOPT.EQ.0) GOTO 110
C
C
C        (D)
C        CHECK REQUESTED EXPONENT VALUES FOR ERRORS
C
      I1 = IEXP+1
      J1 = JEXP+1
      NMIN = MIN0(I1,J1)
      NMAX = MAX0(I1,J1)
      IF (J1.GE.I1) NC = (IEXP+1)*(JEXP+1-IEXP/2)
      IF (J1.LT.I1) NC = (JEXP+1)*(IEXP+1-JEXP/2)
      IF (NC.GT.N.OR.NC.GT.MAXCOF) GOTO 998
      DO 125 K=1,MAXCOF
      IPOWR(K) = 0
  125 JPOWR(K) = 0
C
C        (E)
C        CALL SUBROUTINE SMSRF TO SMOOTH THE DATA Z=F(X,Y)
C
      CALL SMSRF (X,Y,Z,ZNEW,N,IEXP,JEXP,NCOEF,COEF,IPOWR,JPOWR)
      IF (NCOEF.LT.0) GOTO 120
         DO 130 K=1,N
         EPSLON = EPSLON + (Z(K)-ZNEW(K))**2
  130    CONTINUE
         EPSLON = SQRT(EPSLON)/FLOAT(N)
      GOTO 120
  110 DO 100 K=1,N
  100 ZNEW(K) = Z(K)
C
C
C        (F)
C        DETERMINE THE RANGE OF THE Z DATA UNDER CONSIDERATION
C
  120 ZMIN = Z(1)
      ZMAX = ZMIN
         DO 50 K=2,N
         ZMIN = AMIN1(ZMIN,Z(K))
         ZMAX = AMAX1(ZMAX,Z(K))
   50    CONTINUE
      K=0
C
C
C        (G,M)
C        HAS A CONTOUR LIST BEEN GIVEN? . .
C
      K=0
      FN = 1.0
      FK = -1.0
  200 IF (NCNTRS.GT.0) GOTO 180
C
C
C        (H)
C        CALL SUBROUTINE CBVCHK TO VERIFY THAT THE SPECIFIED BASE
C        VALUE IS WITHIN RANGE OF DATA, RESET IF NEEDED
C
      CALL CBVCHK (CLIST(1),CLIST(2),ZMIN,ZMAX,CLNEW)
      IF (CLIST(1).NE.CLNEW) CLIST(1)=CLNEW
C
C        (I)
C        DETERMINE (NEXT) CONTOUR CONSTANT VALUE
C
  210 FK = FK+1.0
      ZCON = FK*FN*CLIST(2) + CLIST(1)
      IF (ZCON.GT.ZMIN.AND.ZCON.LT.ZMAX) GOTO 150
      IF (FN.LT.0.) GOTO 300
      FK = 0.0
      FN = -1.0
      GOTO 210
C
  180 K = K+1
      IF (K.GT.NCNTRS) GOTO 300
      ZCON = CLIST(K)
      IF (ZCON.LT.ZMIN.OR.ZCON.GT.ZMAX) GOTO 200
C
C        (J)
C        CALL SUBROUTINE INTERP TO
C        INTERPOLATE FOR CONTOUR LINE DATA POINTS
C
  150 CALL INTERP (X,Y,ZNEW,N,ZCON,LEDGES,IE,ISMOPT,LAMBDA,
     *               XI,ETA,J,COEF,IPOWR,JPOWR,NCOEF)
C
C        (K,L)
C        ANY DATA POINTS FOUND? . .
C        CALL SUBROUTINE CNTOUR TO SORT THE INTERPOLATED POINTS
C        ON THE CONTOUR LINE AND DRAW IT
C
      IF (J.NE.0) CALL CNTOUR (ZCON,XI,ETA,LAMBDA,J,IBE,ITE)
      GOTO 200
C
  300 RETURN
  997 IERR = 1
      RETURN
  998 IERR = 2
      RETURN
      END
      SUBROUTINE TRIANG (XD,YD,N,L,E,BE,TE)
C
C     ------------------------------------------------------------------
C
C        A SET OF N DATA POINTS ARE KNOWN (X(I),Y(I),I=1,N)  THEY ARE TO
C        BE CONNECTED BY LINES TO FORM A SET OF TRIANGLES (FOR N.LE.
C        MAXPTS).  THE FINAL TRIANGULATION ESTABLISHES A CONVEX POLYGON
C        DEFINED BY LINKED LISTS OF EDGE NUMBERS, END POINTS AND
C        BOUNDARY EDGES.
C
C
C
C     SUBROUTINE INPUT
C       XD   = ARRAY OF ABSCISSAS
C       YD   = ARRAY OF ORDINATES
C       N   = NUMBER OF POINTS IN X AND Y
C
C     SUBROUTINE OUTPUT
C       L   = NUMBER OF EDGES LISTED IN E, BE AND TE
C       E   = LIST OF INDICES OF EACH TRAINGLE EDGE
C       BE  = 1 IF I OF E IS A BOUNDARY EDGE
C       TE  = INDICES OF NEIGHBORING EDGES FOR EACH TRAINGLE
C
C     LOCAL VARIABLES
C       P   = INDICES OF POINTS OUTSIDE THE BOUNDARY
C       J   = NO. OF VALUES IN LIST P
C       B   = INDEX OF POINTS ON THE BOUNDARY ..  INORDER
C       K   = NO. OF POINTS LISTED IN ARRAY B
C       T   = INDICES OF ADJACENT TRIANGLE EDGES
C       M   = NO. OF ROWS USED IN ARRAY T
C       X   = ARRAY OF SCALED X DATA
C       Y   = ARRAY OF SCALED Y DATA
C
C     ------------------------------------------------------------------
C
C
      IMPLICIT INTEGER (P,B)
      INTEGER  T,TE,E
C
      DIMENSION  XD(N),YD(N),X(500),Y(500)
      DIMENSION  P(500),B(500)
      DIMENSION  E(1494,2),BE(1494),TE(1494,4)
      DIMENSION  T(995,3)
C
C     ..DOUBLE PRECISION SPECIFICATION STATEMENTS FOR IBM360
      REAL*8    TERM,DCOMP,D,D1,S,TC
      REAL*8    XP1,X21,YP1,Y21,XP2,X12,YP2,Y12,X1P,Y1P,X2P,Y2P
C
C
C        (A)
C        THE PROCEDURE BEGINS WITH NO BOUNDARY, NO EDGES, AND
C        ALL X-Y DATA POINTS UNDER CONSIDERATION
C        SCALE THE X,Y DATA AND INITIALIZE LOCAL VARIABLES.
C
C
      J = N
      K = 0
      L = 0
      M = 0
      KKNT = 0
      DO 100 JCNT=1,J
  100 P(JCNT) = JCNT
      XMAX = XD(1)
      XMIN = XD(1)
      YMAX = YD(1)
      YMIN = YD(1)
         DO 98 K=2,N
         XMAX = AMAX1(XMAX,XD(K))
         XMIN = AMIN1(XMIN,XD(K))
         YMAX = AMAX1(YMAX,YD(K))
         YMIN = AMIN1(YMIN,YD(K))
   98    CONTINUE
      DLXINV = 1.0/(XMAX-XMIN)
      DLYINV = 1.0/(YMAX-YMIN)
         DO 99 K=1,N
         X(K) = XD(K)*DLXINV
         Y(K) = YD(K)*DLYINV
   99    CONTINUE
C
C
C
C        (B)
C        BEGIN BY TAKING THE LAST PAIR OF POINTS (X,Y(J)) IN THE
C        LIST TO BE THE FIRST BOUNDARY POINT
C
      B(1) = J
      J = J-1
C
C
C
C        (C)
C        FROM THE REAMINING POINTS, FIND THE POINT NEAREST THE FIRST
C
C
      I2 = 1
      I1 = B(1)
      DMIN = (X(I1)-X(1))**2 + (Y(I1)-Y(1))**2
      DO 270 J1=2,J
      DST = (X(I1)-X(J1))**2 + (Y(I1)-Y(J1))**2
      IF (DST.GE.DMIN) GOTO 270
      I2 = J1
      DMIN = DST
  270 C  ONTINUE
C
C        (D)
C        NOW B(1) TO B(I2) IS THE FIRST EDGE.
C        THERE IS ONE EDGE AND TWO BOUNDARY POINTS.
C
      J = J-1
      IF (I2.GT.J) GOTO 275
      DO 274 JCNT=I2,J
      P(JCNT) = P(JCNT+1)
  274 CONTINUE
  275 K = 2
      B(2) = I2
      L = 1
      E(1,1) = MIN0(B(1),B(2))
      E(1,2) = MAX0(B(1),B(2))
C
C
C
C        (E)
C        NOW BEGIN CIRCLING AROUND THE BOUNDARY OF THE POLYGON,
C        CONSIDERING, IN ORDER, EACH BOUNDARY EDGE.  MAINTAIN THE
C        FOLLOWING INDICES -
C          K1 = B ARRAY INDEX OF THE CURRENT EDGE - POINT 1
C          K2 = B ARRAY INDEX OF THE CURRENT EDGE - POINT 2
C          B1,B2 = INDICES OF BOUNDARY POINT COORDINATES
C
      K1 = 0
      KT = 0
   11 K1 = K1+1
      IF (K1.GT.K) K1=1
   12 K2 = K1+1
      IF (K2.GT.K) K2=1
      B1 = B(K1)
      B2 = B(K2)
      KT = KT+1
C
C        (F)
C        CONSIDER THE BOUNDARY EDGE FROM B1 TO B2.  FOR ALL POINTS NOT
C        YET TRIANGULATED (THE J POINTS REMAINING IN P), FIND THE
C        POINT THAT, WHEN TRAINGULATED WITH B1,B2, MINIMIZES THE LENGTH
C        OF THE TWO NEW EDGES TO BE DRAWN.
C
      D1 = 0.
      J1 = 0
      BFLAG = 0
      IF (J.EQ.0) GOTO 6
      DO 1 LJ=1,J
      PJ = P(LJ)
      TERM = (Y(PJ)-Y(B1))*(X(B2)-X(B1))-(X(PJ)-X(B1))*(Y(B2)-Y(B1))
      IF (TERM.LE.0.) GOTO 1
      D = SQRT((X(PJ)-X(B1))**2+(Y(PJ)-Y(B1))**2)                       -
     2   +SQRT((X(PJ)-X(B2))**2+(Y(PJ)-Y(B2))**2)
      IF (J1.NE.0.AND.D1.LT.D) GOTO 1
      J1 = LJ
      D1 = D
    1 CONTINUE
C
C        (G)
C        IF LESS THAN THREE EDGES EXIST (NO TRIANGLE DEFINED YET),
C        THEN THERE ARE NO ADJACENT BOUNDARY POINTS TO BE CONSIDERED.
C        SO GO TO SECTION J.
C
      IF (K.LE.3) GOTO 3
C
C
C        (H)
C        CONSIDER THE ADJACENT BOUNDARY POINT OF THE NEXT EDGE OF THE
C        POLYGON.  CALL ITS INDEX NUMBER K3 AND SEE IF ITS CLOSER TO
C        THE CURRENT EDGE THAN P(J1).
C
    6 K3 = K2+1
      IF (K3.GT.K) K3=1
      PK3 = B(K3)
      TERM = (Y(PK3)-Y(B1))*(X(B2)-X(B1))-(X(PK3)-X(B1))*(Y(B2)-Y(B1))
      IF (TERM.LE.0.) GOTO 2
      D = SQRT((X(PK3)-X(B1))**2+(Y(PK3)-Y(B1))**2)                     -
     2   +SQRT((X(PK3)-X(B2))**2+(Y(PK3)-Y(B2))**2)
      IF (J1.NE.0.AND.D1.LT.D) GOTO 2
      J1 = K3
      D1 = D
      BFLAG = 1
C
C        (I)
C        CONSIDER THE ADJACENT BOUNDARY POINT OF THE PREVIOUS EDGE OF
C        THE POLYGON.  CALL ITS INDEX NUMBER K0 AND SEE IF ITS CLOSER
C        TO THE CURRENT EDGE THAN P(J1) AND B(K3).
C
    2 CONTINUE
      K0 = K1-1
      IF (K0.LT.1) K0=K
      PK0 = B(K0)
      TERM = (Y(PK0)-Y(B1))*(X(B2)-X(B1))-(X(PK0)-X(B1))*(Y(B2)-Y(B1))
      IF (TERM.LE.0.) GOTO 3
      D = SQRT((X(PK0)-X(B1))**2+(Y(PK0)-Y(B1))**2)                     -
     2   +SQRT((X(PK0)-X(B2))**2+(Y(PK0)-Y(B2))**2)
      IF (J1.NE.0.AND.D1.LT.D) GOTO 3
      J1 = K0
      D1 = D
      BFLAG = -1
    3 CONTINUE
C
C        (J)
C        SKIP THE NEXT SECTION IF J1 IS STILL ZERO, SINCE A CANDIDATE
C        POINT FOR TRIANGULATION WITH EDGE B1,B2 WAS NOT FOUND.
C
      IF (J1.EQ.0) GOTO 9
C
C
C
C
C        (K,L)
C        IF THE SEARCH FOR A CANDIDATE POINT HAS ALREADY CONSIDERED EACH
C        BOUNDARY EDGE AT LEAST ONCE (KT.GT.K) OR IF THE BOUNDARY IS
C        BEING CHECKED FOR CONCAVE EDGES (J=0), THEN THE NEXT SECTION
C        (SECTION M) CAN BE OMMITTED.
C
      IF (KT.GT.K.OR.J.EQ.0) GOTO 9
C
C        (M)
C        AT THIS POINT THE USER MAY INSERT ANY ADDITIONAL CONSTRAINT
C        ON THE TRIANGLE TO BE FORMED BY THE POINT PJ1.  IF THE
C        CANDIDATE TRIANGLE FAILS THE TEST, IT IS DELETED FROM
C        CONSIDERATION BY SETTING THE VARIABLE J1 TO ZERO.
C
    9 CONTINUE
C
C
C
C        (N,O)
C        THE NEXT PROCEDURE CHECKS ALL BOUNDARY EDGES OF THE POLYGON
C        FOR INTERSECTION WITH THE CANDIDATE TRIANGLE.  IF ANY EXISTING
C        BOUNDARY EDGE INTERSECTS ANY OF THE EDGES TO BE FORMED BY THE
C        CANDIDATE TRIANGLE, THEN THE CANDIDATE POINT IS REJECTED.  IF
C        BFLAG IS NOT ZERO, THEN THE EDGE DEFINED BY J1=K0 OR J1=K3 IS
C        EXEMPT FORM THIS TEST.
C
C        IF THERE ARE THREE OR LESS EXISTING BOUNDARY EDGES OR IF
C        J1 HAS BEEN SET TO ZERO, THIS TEST IS OMMITTED.
C
      IF (K.LE.3.OR.J1.EQ.0) GOTO 7
      IF (BFLAG.EQ.0) NQ = P(J1)
      IF (BFLAG.EQ.1) NQ = B(K3)
      IF (BFLAG.EQ.-1) NQ = B(K0)
        DO 108 KCNT=1,K
        IF (KCNT.EQ.K1) GOTO 108
        KN = KCNT+1
        IF (KCNT.EQ.K) KN=1
        IF (BFLAG.EQ.-1.AND.(KCNT.EQ.K0.OR.KN.EQ.K0)) GOTO 108
        IF (BFLAG.EQ. 1.AND.(KCNT.EQ.K3.OR.KN.EQ.K3)) GOTO 108
        P1 = B(KCNT)
        P2 = B(KN)
          DO 8 JCNT=1,2
          IF (JCNT.EQ.1.AND.(BFLAG.EQ.0.OR.BFLAG.EQ.1).AND.KCNT.EQ.K0)  -
     *        GOTO 108
          IF (JCNT.EQ.2.AND.(BFLAG.EQ.0.OR.BFLAG.EQ.-1).AND.KCNT.EQ.K2) -
     *        GOTO 108
          BJ = B1
          IF (JCNT.EQ.2) BJ=B2
          XQB = X(NQ)-X(BJ)
          YQB = Y(NQ)-Y(BJ)
          X12 = X(P1)-X(P2)
          Y12 = Y(P1)-Y(P2)
          D = XQB*Y12-YQB*X12
          IF (D.EQ.0.) GOTO 8
          X1B = X(P1)-X(BJ)
          Y1B = Y(P1)-Y(BJ)
          S = (X1B*Y12-Y1B*X12)/D
          IF (S.LT.0..OR.S.GT.1.) GOTO 8
          TC = (XQB*Y1B-YQB*X1B)/D
          IF (TC.LT.0.OR.TC.GT.1.) GOTO 8
          J1 = 0
          GOTO 7
    8     CONTINUE
  108   CONTINUE
    7 CONTINUE
C
C
C        (P,Q)
C        IF J1 IS ZERO, THEN THE CANDIDATE POINT DID NOT PASS THE ABOVE
C        TESTS OR NO POINT WAS FOUND.  IF BFLAG IS NOT ZERO, THEN A
C        POINT ON THE BOUNDARY WAS FOUND.
C
      IF (J1.EQ.0) GOTO 10
      IF (BFLAG) 150,160,4
C
C
C        THE TRIANGULATED POINT IS OUTSIDE THE BOUNDARY.  ESTABLISH TWO
C        NEW EDGES, A NEW BOUNDARY POINT AND DELETE ONE POINT FROM
C        OUTSIDE THE BOUNDARY.
C
C
  160 E(L+1,1) = MIN0(P(J1),B(K1))
      E(L+1,2) = MAX0(P(J1),B(K1))
      E(L+2,1) = MIN0(P(J1),B(K2))
      E(L+2,2) = MAX0(P(J1),B(K2))
      KT = 0
      L = L+2
      M = M+1
      T(M,1) = MIN0(P(J1),B(K1),B(K2))
      T(M,2) = MIDDLE(P(J1),B(K1),B(K2))
      T(M,3) = MAX0(P(J1),B(K1),B(K2))
      IF (K1.EQ.K) GOTO 140
          KM = K
          KP1 = K1+1
  147     B(KM+1) = B(KM)
          KM = KM-1
          IF (KM.GE.KP1) GOTO 147
  140 B(K1+1) = P(J1)
      K = K+1
      J = J-1
      IF (J1.GT.J) GOTO 10
      DO 144 JCNT=J1,J
  144 P(JCNT) = P(JCNT+1)
      GOTO 10
C
C
C        (S)
C        THE TRIANGULATED POINT IS THE NEXT POINT ON THE BOUNDARY.
C        ESTABLISH ONE NEW EDGE (FROM B(K1) TO B(K3)), ONE NEW
C        TRIANGLE (FROM B(K1) TO B(K2) TO B(K3)), AND DELETE ONE POINT
C        FROM THE BOUNDARY (B(K2)).
C
C
    4 E(L+1,1) = MIN0(B(K3),B(K1))
      E(L+1,2) = MAX0(B(K3),B(K1))
      KK = 0
      KKNT = 0
      KT = 0
      L=L+1
      K=K-1
      M = M+1
      T(M,1) = MIN0(B(K1),B(K2),B(K3))
      T(M,2) = MIDDLE(B(K1),B(K2),B(K3))
      T(M,3) = MAX0(B(K1),B(K2),B(K3))
      IF (K2.GT.K) GOTO 155
      DO 151 KCNT=K2,K
  151 B(KCNT) = B(KCNT+1)
  155 IF (K2.EQ.1) K1=K1-1
      GOTO 10
C
C        (R)
C        THE TRIANGULATED POINT IS THE PREVIOUS POINT ON THE BOUNDARY.
C        ESTABLISH A NEW EDGE (FROM B(K0) TO B(K2)), ONE NEW TRIANGLE
C        (FROM B(K0) TO B(K1) TO B(K2)), AND DELETE ONE POINT FROM THE
C        BOUNDARY (B(K1))
C
  150 E(L+1,1) = MIN0(B(K0),B(K2))
      E(L+1,2) = MAX0(B(K0),B(K2))
      KK = 0
      KKNT = 0
      KT = 0
      L = L+1
      K = K-1
      M = M+1
      T(M,1) = MIN0(B(K0),B(K1),B(K2))
      T(M,2) = MIDDLE(B(K0),B(K1),B(K2))
      T(M,3) = MAX0(B(K0),B(K1),B(K2))
      IF (K1.GT.K) GOTO 157
      DO 158 KCNT=K1,K
  158 B(KCNT) = B(KCNT+1)
  157 K1 = K1-1
      IF (K1.LT.1) K1=K
C
C
C        (T)
C        IF THERE ARE ANY POINTS REMAINING OUTSIDE THE BOUNDARY, THEN
C        REPEAT THE PROCEDURE FOR THE NEXT EDGE.
C
C
   10 IF (J.GT.0.AND.J1.NE.0) GOTO 12
      IF (J.GT.0) GOTO 11
C
C
C        (U,V,W)
C        ALL POINTS HAVE BEEN TRIANGULATED.  CHECK THAT ALL BOUNDARY
C        EDGES FORM A CONCAVE POLYGON.
C
      IF (KK.NE.0) GOTO 55
      KK = 1
      KL = 0
   55 KKNT = KKNT+1
      IF (KKNT.GT.N) GOTO 170
    5 KL = KL+1
      K2 = KL+1
      IF (K2.GT.K) K2=1
      K1 = KL-1
      IF (K1.LT.1) K1=K
      PKL = B(KL)
      B1 = B(K1)
      B2 = B(K2)
      TERM = (Y(PKL)-Y(B1))*(X(B2)-X(B1))-(X(PKL)-X(B1))*(Y(B2)-Y(B1))
      IF (TERM.LT.0.) GOTO 11
      IF (KL.LT.K) GOTO 5
C
C
C        (X)
C        THE TRIANGULATION IS COMPLETE AND HAS BEEN CHECKED FOR A
C        CONCAVE BOUNDARY.  NOW IDENTIFY THE BOUNDARY EDGES.
C
C
  170 DO 23 LCNT=1,L
      BE(LCNT) = 0
      KL = 0
   21 KL = KL+1
      IF (E(LCNT,1).NE.B(KL)) GOTO 22
      K1 = KL+1
      IF (K1.GT.K) K1=1
      IF (E(LCNT,2).NE.B(K1)) GOTO 162
      BE(LCNT) = 1
      GOTO 23
  162 K1 = KL-1
      IF (K1.LT.1) K1=K
      IF (E(LCNT,2).NE.B(K1)) GOTO 22
      BE(LCNT) = 1
      GOTO 23
   22 IF (KL.LT.K) GOTO 21
   23 CONTINUE
C
C
C        (Y)
C        FINALLY, ESTABLISH THE INDICES OF ADJACENT EDGES FOR EACH
C        EDGE IN THE TRIANGULATION.  EACH BOUNDARY EDGE WILL HAVE TWO
C        ADJACENT EDGES - EACH INTERIOR EDGE WILL HAVE FOUR.
C
      DO 190 LL  =1,4
      DO 190 LCNT=1,L
  190 TE(LCNT,LL) = 0
      DO 191 MCNT=1,M
         DO 192 LL=1,L
         IF (E(LL,1).EQ.T(MCNT,1).AND.E(LL,2).EQ.T(MCNT,2)) L1=LL
         IF (E(LL,1).EQ.T(MCNT,2).AND.E(LL,2).EQ.T(MCNT,3)) L2=LL
         IF (E(LL,1).EQ.T(MCNT,1).AND.E(LL,2).EQ.T(MCNT,3)) L3=LL
  192    CONTINUE
      LAMBDA = 0
      IF (TE(L1,1).NE.0) LAMBDA=2
      TE(L1,LAMBDA+1) = L2
      TE(L1,LAMBDA+2) = L3
      LAMBDA = 0
      IF (TE(L2,1).NE.0) LAMBDA=2
      TE(L2,LAMBDA+1) = L1
      TE(L2,LAMBDA+2) = L3
      LAMBDA = 0
      IF (TE(L3,1).NE.0) LAMBDA = 2
      TE(L3,LAMBDA+1) = L1
      TE(L3,LAMBDA+2) = L2
  191 CONTINUE
C
      RETURN
      END
      SUBROUTINE SMSRF (X,Y,Z,ZNEW,N,I,J,NCOEF,CNORM,IPOWR,JPOWR)
C
C     ------------------------------------------------------------------
C
C     SUBROUTINE SMSRF PERFORMS THE OPTIONAL SMOOTHING OF DATA BEFORE
C     TRIANGULATION OF THE PLANE IS INITIATED.  THE SURFACE DEFINED BY
C     Z = F(X,Y) IS SMOOTHED VIA A POLYNOMIAL CURVE FIT DEFINED BY A
C     LEAST SQUARES CRITERIA.
C
C
C     ARGUMENTS --
C     (INPUT)
C        X,Y,Z     ARRAYS OF VALUES DEFINING THE KNOWN SURFACE
C                  (POINTS IN SPACE FOR THE FUNCTION Z=F(X,Y))
C        N         THE NUMBER OF POINTS IN X,Y AND Z.
C        I,J       ARE THE EXPONENTS FOR THE SMOOTHING POLYNOMIAL
C                  AS SELECTED BY THE USER.
C     (RETURN)
C        ZNEW      IS THE ARRAY OF SMOOTHED VALUES FOR THE FUNCTION
C                  (ZNEW WILL CONTAIN THE ORIGINAL Z DATA ON RETURN
C                  IF THE SMOOTHING OPERATION FAILS, IN WHICH CASE
C                  NCOEF WILL BE SET TO -1).
C        NCOEF     IS THE NUMBER OF TERMS IN THE POLYNOMIAL RESULTING
C                  FROM THE VALUES OF I AND J.  NCOEF MUST BE LESS THAN
C                  OR EQUAL TO BOTH N AND MAXCOF.
C        C         IS THE ARRAY OF NCOEF COMPUTED COEFFICIENTS
C        IPOWR     THE ARRAY OF I EXPONENTS FOR EACH TERM
C        JPOWR     THE ARRAY OF J EXPONENTS FOR EACH TERM
C                  (EACH ELEMENT OF C, IPOWR AND JPOWR IS ASSOCIATED
C                  WITH THE NCOEF TERMS OF THE POLYNOMIAL, IN ORDER)
C
C     ------------------------------------------------------------------
C
C
C
      DIMENSION  X(N),Y(N),Z(N),ZNEW(N)
      DIMENSION  IPOWR(23),JPOWR(23),C(23),CNORM(23),AVE(23)
      DIMENSION  IP(23),XX(23),H(23)
      DIMENSION  B(500),AM(500,23)
C
      DATA       IA /500/
C
C
C        (A)
C        INITIALIZE LOCAL VARIABLES AND RANGE CHECK
C
      REALN = FLOAT(N)
      IF(I.LT.1) I = 1
      IF(J.LT.1) J = 1
      I1 = I+1
      J1 = J+1
      NCOEF = 0
C
C
C        (B)
C        DETERMINE THE X AND Y EXPONENTS TO BE USED
C        SAVE THEM IN ARRAYS IPOWR AND JPOWR
C
      NCOEF = 0
      K = MAX0(I1,J1)
      IF (K.EQ.0) GOTO 950
          DO 180 II=1,I1
          KI1 = K-II+1
          L = MIN0(KI1,J1)
             DO 181 JJ=1,L
             NCOEF = NCOEF+1
             IPOWR(NCOEF) = II-1
             JPOWR(NCOEF) = JJ-1
  181        CONTINUE
  180     CONTINUE
C
C
C        (C)
C        USING THE EXPONENT LISTS FROM ABOVE AND THE
C        KNOWN XY DATA POINTS, CONSTRUCT THE MATRIX AM
C
      DO 182 KCOL=1,NCOEF
      IEX = IPOWR(KCOL)
      JEX = JPOWR(KCOL)
          DO 284 KROW=1,N
          X2 = X(KROW)
          IF (X2.EQ.0.0) X2=1.0
          XP = X2**IEX
          Y2 = Y(KROW)
          IF (Y2.EQ.0.0) Y2=1.0
          YP = Y2**JEX
              AM(KROW,KCOL) = XP*YP
  284     CONTINUE
  182 CONTINUE
      KROW = NCOEF
C
C
C        (D)
C        NORMALIZE EACH VALUE IN EACH COLUMN OF AM BY THE COLUMN AVERAGE
C
      AVE(1) = 1.0
      DO 403 L1 = 2,NCOEF
      AVE(L1) = 0.0
      DO 402 L2 = 1,N
  402 AVE(L1) = AVE(L1) +  ABS(AM(L2,L1))
      AVE(L1) = AVE(L1)/REALN
      IF (AVE(L1) .EQ. 0.)  AVE(L1) = 1.0
      DO 404 L2 = 1,N
  404 AM(L2,L1) = AM(L2,L1)/AVE(L1)
  403 CONTINUE
C
C
C
C        (E,F,G)
C        USE IMSL ROUTINE LLSQF TO SOLVE (VIA LEAST-SQUARES)
C        THE SYSTEM  AM*C = Z  FOR MATRIX C
C
C
      M = N
      IER = 0
      KBASIS = NCOEF
      TOL = 0.0
      DO 222 KK=1,N
      B(KK) = Z(KK)
  222 CONTINUE
      CALL LLSQF (AM,IA,M,NCOEF,B,TOL,KBASIS,XX,H,IP,IER)
      IF (IER.NE.0) GOTO 950
C
C
C        (H)
C        DIVIDE OUT THE SCALE FACTOR FROM THE SOLUTION
C        MATRIX AND ESTABLISH THE COEFFICIENTS
C
      DO 905 L3 = 1,NCOEF
      C(L3) = XX(L3)
      CNORM(L3) = C(L3)/AVE(L3)
  905 CONTINUE
C
C
C        (I)
C        ESTABLISH THE NEW Z VALUES BY
C        EVALUATING THE POLYNOMIAL FOR EACH KNOWN X-Y PAIR
C
      DO 934 L3=1,N
      ZNEW(L3) = -1.0*POLYX2(0.0,X(L3),Y(L3),CNORM,IPOWR,JPOWR,NCOEF)
  934 CONTINUE
      RETURN
C
C
C
C        (J)
C        ERROR RETURN,  SET NCOEF TO -1 AND
C        SEND BACK OLD Z VALUES TO CALLING PROGRAM
C
  950 DO 960 L1=1,N
  960 ZNEW(L1) = Z(L1)
      NCOEF = -1
      RETURN
C
C
      END
      FUNCTION MIDDLE(I,J,K)
C
C
C        THIS FUNCTION SUBPROGRAM IS USED BY THE TRIANGULATION ALGORITHM
C        TO FIND THE MIDDLE VALUE OF THE THREE INTEGER ARGUMENTS (THE
C        VALUE WHICH IS NEITHER A MINIMUM OR A MAXIMUM).  I, J AND K ARE
C        ARE ASSUMED TO BE DISCRETE VALUES WITH NO TWO EQUAL.
C
C
      IF (J.LT.I.AND.I.LT.K) GOTO 100
      IF (K.LT.I.AND.I.LT.J) GOTO 100
      IF (I.LT.J.AND.J.LT.K) GOTO 200
      IF (K.LT.J.AND.J.LT.I) GOTO 200
      MIDDLE = K
      RETURN
  100 MIDDLE = I
      RETURN
  200 MIDDLE = J
      RETURN
      END
      FUNCTION POLYX2 (Z,X,Y,C,IPOWR,JPOWR,NCOEF)
C
C
C        POLYX2 IS THE POLYNOMIAL EVALUATION FUNCTION USED WHEN THE
C        SMOOTHING OPTION HAS BEEN INVOKED.  X AND Y LISTS ARE THE
C        KNOWN VALUES OF THE INDEPENDENT VARIABLES.  C IS THE LIST OF
C        COEFFICIENTS FOR EACH TERM.  IPOWR AND JPOWR ARE THE EXPONENTS
C        FOR EACH TERM AND N IS THE NUMBER OF TERMS IN THE POLYNOMIAL.
C        Z IS AN OFFSET TERM WHEN EVALUATING FOR A CONSTANT X VALUE.
C
C
      DIMENSION  IPOWR(23),JPOWR(23),C(23)
C
      IF ((X.EQ.0.0) .OR. (Y.EQ.0.0)) GOTO 140
         POLYX2 = C(1)
         IF (NCOEF.EQ.1) GOTO 130
            DO 120 II=2,NCOEF
            POLYX2 = POLYX2 + (X**IPOWR(II))*(Y**JPOWR(II))*C(II)
  120       CONTINUE
  130    POLYX2 = Z-POLYX2
      RETURN
C
C
  140 POLYX2 = C(1)
      IF (NCOEF.EQ.1) GOTO 160
         DO 150 II=2,NCOEF
         TERM1 = 1.0
         IF (IPOWR(II).NE.0) TERM1 = X**IPOWR(II)
         TERM2 = 1.0
         IF (JPOWR(II).NE.0) TERM2 = Y**JPOWR(II)
         POLYX2 = POLYX2 + TERM1*TERM2*C(II)
  150    CONTINUE
  160 POLYX2 = Z-POLYX2
      RETURN
      END
      SUBROUTINE CBVCHK (ZZERO,DELZ,ZMIN,ZMAX,ZZNEW)
C
C     ------------------------------------------------------------------
C        CONTOUR BASE VALUE CHECKING ROUTINE
C        *       *    *     **  *
C
C     THIS SUBROUTINE SHIFTS THE BASE VALUE (ZZERO) UNTIL IT FALLS
C     WITHIN THE RANGE OF DATA FOR THIS CONTOUR (I.E. BETWEEN ZMIN
C     AND ZMAX).  THE SHIFTED VALUE (THE NEW STARTING BASE VALUE) IS
C     RETURNED TO CALLER AS ZZNEW.  THE USER SHIFT INCREMENT COMES
C     INTO CBVCHK AS DELZ FOR Z CONTOURS.
C
C     ARGUMENTS -
C        ZZERO     = BASE VALUE (INPUT)
C        DELZ      = INCREMENT VALUE (INPUT)
C        ZMIN,ZMAX = RANGE OF Z DATA (INPUT)
C        ZZNEW     = NEW BASE VALUE, MAY OR MAY NOT BE
C                    THE SAME AS ZZERO (RETURN)
C
C
C     ------------------------------------------------------------------
C
C
      IF (ZMIN.EQ.ZMAX) GOTO 999
      ZZNEW = ZZERO
      IF (ZMIN.LE.ZZNEW.AND.ZZNEW.LE.ZMAX) GOTO 999
    2 IF (ZZNEW.GE.ZMAX) GOTO 1
      ZZNEW = ZZNEW + DELZ
      IF (ZMIN.LE.ZZNEW.AND.ZZNEW.LE.ZMAX) GOTO 999
      GOTO 2
C
    1 ZZNEW = ZZERO
    4 IF (ZZNEW.LE.ZMIN) GOTO 999
      ZZNEW = ZZNEW - DELZ
      IF (ZMIN.LE.ZZNEW.AND.ZZNEW.LE.ZMAX) GOTO 999
      GOTO 4
C
  999 RETURN
      END
      SUBROUTINE INTERP (X,Y,U,N,ZCON,LEDGES,IE,ISMOPT,LAMBDA,XI,
     *                   ETA,J,C,IPOWR,JPOWR,NCOEF)
C     ------------------------------------------------------------------
C
C     SUBROUTINE INTERP IS GIVEN A CONSTANT U VALUE (BIGU) FOR WHICH
C     THE CONTOUR LINE IS TO BE DRAWN.  CHECK ALL GIVEN TRIANGLE EDGES,
C     (ARRAY IE) AND CHECK THE VALUES OF U AT THE ENDPOINTS.
C     INTERPOLATE FOR ALL POSSIBLE VALUES ON THE TRIANGLE EDGES.
C     IF ISMOPT = 0, THEN USE A LINEAR INTERPOLATION, IF ISMOPT NOT ZERO
C     THEN EVALUATE FOR A NON-LINEAR SURFACE USING THE COEFFICIENTS
C     FROM SMSRF AND FUNCTION SUBROUTINE POLYX.
C
C        X,Y,U     = DEPENDENT AND INDEPENDENT VALUES FOR
C                    THE RELATION U=F(X,Y)  (INPUT)
C        ZCON      = CONSTANT VALUE OF Z FOR WHICH INTERPOLATION
C                    IS REQUIRED  (INPUT)
C        LEDGES    = NO. OF EDGES IN THE TRIANGULATION  (INPUT)
C        IE        = EDGE ENDPOINT INDICES FROM TRIANGULATION (INPUT)
C        ISMOPT    = SMOOTHING OPTION FLAG,  0=OFF, 1=ON,  (INPUT)
C        LAMBDA    = INDEX OF EDGES FOR INTERPOLATED POINTS  (RETURN)
C        XI        = LIST OF X-COORDINATES OF INTERPOLATED POINTS
C        ETA       = LIST OF Y-COORDINATES OF INTERPOLATED POINTS
C        J         = NUMBER OF VALUES IN XI, ETA LISTS
C                    (XI, ETA AND J ARE RETURNED)
C        C         = LIST OF COEFFICIENTS OF EACH TERM OF THE EQUATION.
C        IPOWR,JPOWR ARE THE LIST OF EXPONENTS FOR EACH TERM OF
C        THE POLYNOMIAL USED TO SMOOTH THE DATA  (INPUT).
C        NCOEF     = NUMBER OF TERMS IN THE POLYNOMIAL
C                    (IPOWR,JPOWR,C, AND NCOEF ARE INPUT)
C
C
C     ------------------------------------------------------------------
C
C
      DIMENSION  X(N),Y(N),U(N)
      DIMENSION  IE(1494,2),XI(1494),ETA(1494),LAMBDA(1494)
      DIMENSION  IPOWR(23),JPOWR(23),C(23)
C
C
      IF (NCOEF.LT.1) ISMOPT=0
      J = 0
C
      DO 1 LCNT=1,LEDGES
C
C        (A)
C        DETERMINE X,Y,Z FOR THE ENDPOINTS OF THE NEXTEDGE - ORDER THEM
C
      I1 = IE(LCNT,1)
      I2 = IE(LCNT,2)
      X1 = X(I1)
      X2 = X(I2)
      Y1 = Y(I1)
      Y2 = Y(I2)
      U1 = U(I1)
      U2 = U(I2)
C
C        (B)
C        FUNCTION VALUES EQUAL AT ENDPOINTS OR
C        CONSTANT ZC NOT BETWEEN THEM? . .
C
      IF (U1.EQ.U2) GOTO 1
      IF (U1.LT.U2) GOTO 100
      TEMP = U2
      U2 = U1
      U1 = TEMP
      TEMP = X2
      X2 = X1
      X1 = TEMP
      TEMP = Y2
      Y2 = Y1
      Y1 = TEMP
  100 IF (ZCON.LT.U1.OR.U2.LT.ZCON) GOTO 1
      IF (U2.EQ.ZCON) U2 = 1.000001 * ZCON
      J = J+1
C
C        (C)
C        HAS DATA BEEN SMOOTHED? . .
C        IF NOT, GOTO SECTION E (STATEMENT LABEL 101)
C
      IF (ISMOPT.EQ.0) GOTO 101
C
C        (D,F)
C        NON-LINEAR INTERPOLATION IS REQUIRED
C        ON THIS EDGE OVER THE Z-SURFACE
C
      F1 = POLYX2 (ZCON,X1,Y1,C,IPOWR,JPOWR,NCOEF)
C
      DO 220 K=1,10
      XN = (X1+X2)*0.5
      YN = (Y1+Y2)*0.5
      FN = POLYX2 (ZCON,XN,YN,C,IPOWR,JPOWR,NCOEF)
      IF (FN.EQ.0.) GOTO 132
      IF (FN.LT.0..AND.F1.LT.0.) GOTO 235
      IF (FN.GT.0..AND.F1.GT.0.) GOTO 235
      X2 = XN
      Y2 = YN
      GOTO 220
  235 X1 = XN
      Y1 = YN
  220 CONTINUE
  132 XI(J)  = (X1+X2)*0.5
      ETA(J) = (Y1+Y2)*0.5
      GOTO 200
C
C
C        (E,F)
C        LINEAR INTERPOLATION IS REQUIRED
C        FOR THIS EDGE OVER THE Z-SURFACE
C
  101 T1 = (U2-ZCON)/(U2-U1)
      T2 = (ZCON-U1)/(U2-U1)
      XI(J) = T1*X1+T2*X2
      ETA(J)= T1*Y1+T2*Y2
  200 LAMBDA(J) = LCNT
    1 CONTINUE
      RETURN
      END
      SUBROUTINE CNTOUR (ZCON,XI,ETA,LAMBDA,J,IBE,ITE)
C
C     ------------------------------------------------------------------
C
C     A SET OF J INTERPOLATED POINTS FOR Z=ZCON  (XI(I),ETA(I) ON EDGE
C     LAMBDA(I) FOR I=1,J), THE CONTOUR LINES MUST NOW BE DRAWN.  THERE
C     MAY BE SEVERAL LINES, EITHER OPEN OR CLOSED CONTOURS.  THIS
C     ALGORITHM WILL USE THE TRIANGULATION RELATIONSHIPS TO SORT OUT
C     EACH LINE IN ORDER.  AS EACH CONTOUR LINE IS ESTABLISHED, USER
C     SUPPLIED PROGRAM CNTCRV IS CALLED TO OUTPUT IT TO THE GRAPHICS
C     DEVICE BEING USED.
C
C
C
C     ARGUMENTS (ALL ARE INPUTS) -
C        ZCON    = CONSTANT VALUE OF Z UNDER CONSIDERATION
C        XI(J)   = ARRAY OF X COORDIANTES OF INTERPOLATED POINTS
C        ETA(J)  = ARRAY OF Y COORDIANTES OF INTERPOLATED POINTS
C        LAMBDA(J) = ARRAY OF EDGE NUMBERS FOR J-TH INTERPOLATED POINT
C        J       = NUMBER OF POINTS IN THE LIST OF INTERPOLATED POINTS
C        IBE     = THE LIST OF BOUNDARY EDGES TAKEN FROM THE TRIANGULATION
C        ITE     = LINKED LIST OF INDICES OF ADJACENT EDGES PROVIDED
C                  BY THE TRIANGULATION PROCEDURE.
C
C
C     ------------------------------------------------------------------
C
C
      DIMENSION  XI(1494),ETA(1494),LAMBDA(1494),IBE(1494),XX(1494),
     *           YY(1494),ITE(1494,4)
C
C
C
C        (A)
C        INITIALIZE LOCAL VARIABLES
C
      IF (J.EQ.0) RETURN
   10 J1 = 0
C
C        (B,C)
C        SEARCH THE LIST OF EDGES FOR A BOUNDARY EDGE (BE(I)=1)
C
    1 J1 = J1+1
      L1 = LAMBDA(J1)
      IF (IBE(L1).EQ.1) GOTO 2
      IF (J1.LT.J) GOTO 1
      GOTO 11
C     SEARCH FOR A BOUNDARY EDGE AND PUT IT AT THE TOP OF THE LIST.
C
C        (D)
C        PUT THIS INTERPOLATED POINT AT THE TOP OF THE
C        LIST FOR THIS CONTOUR,  SET J1
C
    2 IF (J1.EQ.J) GOTO 3
      XI(J+1) = XI(J1)
      ETA(J+1) = ETA(J1)
      LAMBDA(J+1) = LAMBDA(J1)
      DO 101 JCNT = J1,J
      XI(JCNT) = XI(JCNT+1)
      ETA(JCNT) = ETA(JCNT+1)
  101 LAMBDA(JCNT) = LAMBDA(JCNT+1)
C
C        (E)
C        SEARCH THE REMAINING POINTS FOR AN ADJACENT (COMMON) EDGE
C
    3 J1BIG = J
      LCNT = L1
    6 J1BIG = J1BIG-1
      J1 = 0
    5 J1 = J1+1
      L1 = LAMBDA(J1)
      DO 102 I=1,4
      IF (L1.EQ.ITE(LCNT,I)) GOTO 4
  102 CONTINUE
C        (F)
C        ERROR - THERE IS NO NEXT POINT.
      IF (J1.LT.J1BIG) GOTO 5
      GOTO 800
C
C        (G)
C        PUT THIS POINT AT THE TOP OF THE
C        LIST.  CONTINUE IF ITS NOT A BOUNDARY EDGE.
C
    4 XI(J+1) = XI(J1)
      ETA(J+1) = ETA(J1)
      LAMBDA(J+1) = LAMBDA(J1)
      DO 103 JCNT = J1,J
      XI(JCNT) = XI(JCNT+1)
      ETA(JCNT) = ETA(JCNT+1)
  103 LAMBDA(JCNT) = LAMBDA(JCNT+1)
      LCNT = L1
      IF (IBE(L1).NE.1) GOTO 6
C
C        (H)
C        DRAW THE OPEN CONTOUR LINE THROUGH THE POINTS
C        XI(J1),ETA(J1) ...... XI(J1+1),ETA(J1+1) ...... XI(J),ETA(J)
C        THEN RESET J AND CONTINUE
C
C
C     ----------------------------------------------------------------
      NPOINT = J-J1BIG+1
      IF (NPOINT.LE.1) GOTO 300
      CALL CNTCRV (XI(J1BIG),ETA(J1BIG),NPOINT,ZCON)
C     ----------------------------------------------------------------
C
  300 J=J1BIG - 1
C
C        (I)
C        ARE THERE ANY MORE POINTS LEFT? . .
C
      IF (J) 800,800,10
C
C
C        (J)
C        NOW DRAW INTERNAL LINES (CLOSED CONTOURS THAT DO NOT START
C        OR STOP AT BOUNDARY EDGES).  THE POINT AT J1BIG=J IN
C        THE LIST IS CHOSEN TO START THE CONTOUR.
C
   11 J1BIG = J+1
      LCNT = LAMBDA(J)
C
C        (K,M,P)
C        FIND THE NEXT POINT FOR THIS CONTOUR (ON AN EDGE WITH A COMMON
C        END POINT).  PUT IT AT THE TOP OF THE LIST, AND REPEAT UNTIL
C        NO MORE COMMON EDGES REMAIN FOR THIS LINE.
C
   16 J1BIG = J1BIG-1
      J1 = 0
      IF (J1BIG.GT.J) J1=1
   15 J1 = J1+1
      L1 = LAMBDA(J1)
      DO 104 I=1,4
      IF (L1.EQ.ITE(LCNT,I)) GOTO 14
  104 CONTINUE
      IF (J1.LT.J1BIG) GOTO 15
C        (L)
C        OTHERWISE, NO ADJACENT EDGE WAS FOUND.
C        THIS CONTOUR LINE IS COMPLETE, GO DRAW IT.
      GOTO 17
C
   14 XI(J+1) = XI(J1)
      ETA(J+1) = ETA(J1)
      LAMBDA(J+1) = LAMBDA(J1)
      DO 105 JCNT = J1,J
      XI(JCNT) = XI(JCNT+1)
      ETA(JCNT) = ETA(JCNT+1)
  105 LAMBDA(JCNT) = LAMBDA(JCNT+1)
      LCNT = L1
      IF (J1BIG.NE.1) GOTO 16
C
C
C
C        (O)
C        DRAW THE CLOSED CONTOUR LINE, THE INTERPOLATED LINE THROUGH
C        XI(J1),ETA(J1) ...... XI(J),ETA(J) ...... XI(J1),ETA(J1)
C
C
C
C     -----------------------------------------------------------------
   17 JJ = J1BIG
      IF (J1BIG.NE.1) JJ = J1BIG+1
      KNT = 0
      DO 510 KK = JJ,J
      KNT = KNT+1
      XX(KNT) = XI(KK)
      YY(KNT) = ETA(KK)
  510 CONTINUE
      XX(KNT+1) = XX(1)
      YY(KNT+1) = YY(1)
      NPOINT = KNT+1
      CALL CNTCRV (XX(1),YY(1),NPOINT,ZCON)
C     -----------------------------------------------------------------
C
C
C        (P)
C        RESET J.  ESTABLISH THE NEXT CONTOUR LINE FOR REMAINING POINTS
C        OR QUIT THE PROCEDURE IF NO MORE POINTS REMAIN.
C
      J = J1BIG - 1
      IF (J.GT.0) J=J+1
      IF (J) 800,800,11
  800 RETURN
      END
